Code implemented in R to generate the figures of the Simonin et al. paper on the wheat microbiome

Figure 1

Load the diversity and structure data for the 8 wheat genotypes in one soil (FR2) - For detailed script on how to generate the diversity indices and NMDS, see Figure 3

Div <- read.table("Diversity&Structure_data_8genotypes_FR2soil.txt", header=TRUE, check.names = FALSE)

Graph ESV Richness - 8 genotypes (Figure 1A)

library(Rmisc)
library(ggplot2)
color= c("#4CB5F5",  "#D70026", "#EDB83D",  "navy",  "#F77604", "#B3C100", "#EC96A4", "#5BC8AC")
Total <- summarySE(Div, measurevar="S", groupvars=c("Genotype"), na.rm = TRUE)
Bac <- summarySE(Div, measurevar="S_16S", groupvars=c("Genotype"), na.rm = TRUE)
Euk <- summarySE(Div, measurevar="S_18S", groupvars=c("Genotype"), na.rm = TRUE)
p2=ggplot() + geom_point(data=Total, aes(x=reorder(Genotype, S), y=S, color=Genotype), size=5, shape=15) + geom_point(data=Bac, aes(x=Genotype, y=S_16S, color=Genotype), size=4, shape=5) + geom_point(data=Euk, aes(x=Genotype, y=S_18S, color=Genotype), size=4, shape=1) + geom_errorbar(data=Total, aes(x=Genotype, ymin=S-se, ymax=S+se, color=Genotype), width=.2)+ geom_errorbar(data=Bac, aes(x=Genotype, ymin=S_16S-se, ymax=S_16S+se, color=Genotype), width=.2) + geom_errorbar(data=Euk, aes(x=Genotype, ymin=S_18S-se, ymax=S_18S+se, color=Genotype), width=.2)+ theme_classic()+xlab("Genotypes")+ylab("Observed ESVs Richness")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

Graph Phylogenetic diversity - 8 genotypes (Figure 1B)

color= c("#4CB5F5",  "#D70026", "#EDB83D",  "navy",  "#F77604", "#B3C100", "#EC96A4", "#5BC8AC")
Total <- summarySE(Div, measurevar="PD", groupvars=c("Genotype"), na.rm = TRUE)
Bac <- summarySE(Div, measurevar="PD_16S", groupvars=c("Genotype"), na.rm = TRUE)
Euk <- summarySE(Div, measurevar="PD_18S", groupvars=c("Genotype"), na.rm = TRUE)
p2=ggplot() + geom_point(data=Total, aes(x=reorder(Genotype, PD), y=PD, color=Genotype), size=5, shape=15) + geom_point(data=Bac, aes(x=Genotype, y=PD_16S, color=Genotype), size=4, shape=5) + geom_point(data=Euk, aes(x=Genotype, y=PD_18S, color=Genotype), size=4, shape=1) + geom_errorbar(data=Total, aes(x=Genotype, ymin=PD-se, ymax=PD+se, color=Genotype), width=.2)+ geom_errorbar(data=Bac, aes(x=Genotype, ymin=PD_16S-se, ymax=PD_16S+se, color=Genotype), width=.2) + geom_errorbar(data=Euk, aes(x=Genotype, ymin=PD_18S-se, ymax=PD_18S+se, color=Genotype), width=.2)+ theme_classic()+xlab("Genotypes")+ylab("Faith's Phylogenetic Diversity")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

Graph correlation between Prokaryotic and Eukaryotic richness - 8 genotypes (Figure 1C)

color= c("#4CB5F5",  "#D70026", "#EDB83D",  "navy",  "#F77604", "#B3C100", "#EC96A4", "#5BC8AC")
p2=ggplot(data=Div) + geom_point(aes(x=S_16S, y=S_18S, color=Genotype), alpha=0.9, size=4) + theme_classic()+xlab("Prokaryotes ESV Richness")+ylab("Eukaryotes  ESV Richness")+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+scale_color_manual(values=color)+geom_smooth(aes(x=S_16S, y=S_18S),method = "lm", se=FALSE, color="black")+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

NMDS genotype effect on total microbiome structure - 8 genotypes (Figure 1D)

color= c("#4CB5F5",  "#D70026", "#EDB83D",  "navy",  "#F77604", "#B3C100", "#EC96A4", "#5BC8AC")
p1=ggplot(data=Div, aes(x=NMDS1, y=NMDS2,color=Genotype))+geom_point(shape=5, size=3.5)+theme_classic(base_size = 15)+xlab("NMDS1")+ylab("NMDS2")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p1

Figure 2

Load the diversity and structure data for the 8 soils (3 wheat genotypes studied) - For detailed script on how to generate the diversity indices and NMDS, see Figure 3

Div <- read.table("Diversity&Structure_data_8soils.txt", header=TRUE, check.names = FALSE)

Graph ESV Richness - 8 soils (Figure 2A)

library(Rmisc)
library(ggplot2)
color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
Total <- summarySE(Div, measurevar="S", groupvars=c("Soils"), na.rm = TRUE)
Bac <- summarySE(Div, measurevar="S_16S", groupvars=c("Soils"), na.rm = TRUE)
Euk <- summarySE(Div, measurevar="S_18S", groupvars=c("Soils"), na.rm = TRUE)
p2=ggplot() + geom_point(data=Total, aes(x=reorder(Soils, S), y=S, color=Soils), size=5, shape=15) + geom_point(data=Bac, aes(x=Soils, y=S_16S, color=Soils), size=4, shape=5) + geom_point(data=Euk, aes(x=Soils, y=S_18S, color=Soils), size=4, shape=1) + geom_errorbar(data=Total, aes(x=Soils, ymin=S-se, ymax=S+se, color=Soils), width=.2)+ geom_errorbar(data=Bac, aes(x=Soils, ymin=S_16S-se, ymax=S_16S+se, color=Soils), width=.2) + geom_errorbar(data=Euk, aes(x=Soils, ymin=S_18S-se, ymax=S_18S+se, color=Soils), width=.2)+ theme_classic()+xlab("Soils")+ylab("Observed ESVs Richness")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

Graph phylogenetic diversity - 8 soils (Figure 2B)

library(Rmisc)
library(ggplot2)
color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
Total <- summarySE(Div, measurevar="PD", groupvars=c("Soils"), na.rm = TRUE)
Bac <- summarySE(Div, measurevar="PD_16S", groupvars=c("Soils"), na.rm = TRUE)
Euk <- summarySE(Div, measurevar="PD_18S", groupvars=c("Soils"), na.rm = TRUE)
p2=ggplot() + geom_point(data=Total, aes(x=reorder(Soils, PD), y=PD, color=Soils), size=5, shape=15) + geom_point(data=Bac, aes(x=Soils, y=PD_16S, color=Soils), size=4, shape=5) + geom_point(data=Euk, aes(x=Soils, y=PD_18S, color=Soils), size=4, shape=1) + geom_errorbar(data=Total, aes(x=Soils, ymin=PD-se, ymax=PD+se, color=Soils), width=.2)+ geom_errorbar(data=Bac, aes(x=Soils, ymin=PD_16S-se, ymax=PD_16S+se, color=Soils), width=.2) + geom_errorbar(data=Euk, aes(x=Soils, ymin=PD_18S-se, ymax=PD_18S+se, color=Soils), width=.2)+ theme_classic()+xlab("Soils")+ylab("Faith's Phylogenetic Diversity")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

Graph correlation between Prokaryotic and Eukaryotic richness - 8 soils (Figure 2C)

color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
p2=ggplot(data=Div) + geom_point(aes(x=S_16S, y=S_18S, color=Soils), alpha=0.9, size=4) + theme_classic()+xlab("Prokaryotes ESV Richness")+ylab("Eukaryotes ESV Richness")+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+scale_color_manual(values=color)+geom_smooth(aes(x=S_16S, y=S_18S),method = "lm", se=FALSE, color="black")+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

NMDS soil effect on total microbiome structure - 8 soils (Figure 2D)

library(ggplot2)
color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
p1=ggplot(data=Div, aes(x=NMDS1, y=NMDS2,color=Soils, shape=Genotype))+geom_point(size=3)+theme_classic(base_size = 15)+xlab("NMDS1")+ylab("NMDS2")+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+scale_color_manual(values=color)+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p1

Figure 3

Load merged 16S and 18S SV table

SV_use1 <- read.table("ESV_table_16S&18Smerged.txt", header=TRUE, check.names = FALSE, sep="\t")

Subset just european soils

Sols=subset(SV_use1, Huit_sols=="8_sols")
SV_use1=subset(Sols, Type=="DNA")
SV_use1=subset(SV_use1, Continent=="Europe")
dim(SV_use1)
## [1]   28 6278

Create filtered matrix for NMDS and diversity indices calculations

matrix<-SV_use1[c(14:6278)]
dim(matrix)
## [1]   28 6265
matrix_use<-matrix[,colSums(matrix)>=1]
dim(matrix_use)
## [1]   28 2552

Ordination: NMDS with Bray-Curtis distances

library(vegan)
NMDS <- metaMDS(matrix_use, distance = "bray", trymax = 100)
NMDS
stressplot(NMDS)

##Extract scores for sites (samples)
NMDSsites=scores(NMDS, display="sites")
SV_use1=cbind(SV_use1,NMDSsites)

Graph effect agricultural pratices in European soils on microbiome structure - NMDS (Figure 3B)

library(ggplot2)
color= c("#4CB5F5", "#D70026")
p1=ggplot(data=SV_use1, aes(x=NMDS1, y=NMDS2,color=Practices, shape=Soils))+geom_point(size=3)+theme_classic(base_size = 15)+xlab("NMDS1")+ylab("NMDS2")+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+scale_color_manual(values=color)+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p1

Calculate Diversity index (S, Shannon, Simpson, Evenness)

Shannon <- diversity(matrix_use)
simp <- diversity(matrix_use, "simpson")
# Species richness (S) and Pielou's evenness (J):
S <- specnumber(matrix_use) 
J <- Shannon/log(S)

divtable= cbind(S,Shannon,J,simp)
SV_use1=cbind(SV_use1,divtable)

Graph effect agricultural pratices in European soils on microbiome ESV richness (Figure 3A)

library(Rmisc)
library(ggplot2)
Diversity_stat <- summarySE(SV_use1, measurevar="S", groupvars=c("Type", "Practices", "Soils"), na.rm = TRUE)
color= c("#4CB5F5", "#D70026")
p2=ggplot() + geom_point(data=Diversity_stat, aes(x=Soils, y=S, color=Practices), size=4)+geom_point(data=SV_use1, aes(x=Soils, y=S, color=Practices), alpha=0.8) + geom_errorbar(data=Diversity_stat, aes(x=Soils, ymin=S-se, ymax=S+se, color=Practices), width=.2) + theme_classic()+xlab("Soils")+ylab("Observed ESVs Richness")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

Figure 4 - Bubble plot “Most abundant clades”

Load the relative abundance data of the most abundant clades (Prokaryotes and Eukaryotes)

Fam <- read.table("Family_bubbleplot.txt", header=TRUE, check.names = FALSE, sep="\t")

Bubble plot of most abundant clades (Figure 4)

library(ggplot2)
mycolors <- scale_color_manual(values = c("black","pink", "orange", "yellow", "purple"))
p2=ggplot(data=Fam) + geom_point(aes(x=Rel_abund, y=(reorder (Clade, Rel_abund)), size=number_SV, color=Group), alpha=0.8) + theme_classic()+xlab("Relative Abundance")+ylab("Clades")+ theme(axis.title = element_text(color="black", size=10, face="bold"))+ theme(axis.text = element_text(color="black", size=10, face="bold"))+ theme(legend.text = element_text(colour="black", size = 10, face = "bold"))+ theme(legend.title = element_text(colour="black", size=10, face="bold"))+
    scale_size_continuous(name = "Number of ESVs", breaks=c(25, 50, 100, 200, 300),
                          limits = c(4, 300),
                          range = c(0, 10) )+ scale_x_continuous(labels = scales::percent_format(accuracy = 1))+mycolors
p2

Figure 5

Load the diversity and structure data for the 8 soils (3 wheat genotypes studied)

Div <- read.table("Diversity&Structure_data_8soils.txt", header=TRUE, check.names = FALSE)

Number of core taxa per sample (Figure 5C)

library(Rmisc)
library(ggplot2)
Diversity_stat <- summarySE(Div, measurevar="core_count", groupvars=c( "Soils"), na.rm = TRUE)
color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
p2=ggplot() + geom_point(data=Diversity_stat, aes(x=reorder(Soils, core_count),, y=core_count, color=Soils), size=4)+geom_point(data=Div, aes(x=Soils, y=core_count, color=Soils), alpha=0.8) + geom_errorbar(data=Diversity_stat, aes(x=Soils, ymin=core_count-se, ymax=core_count+se, color=Soils), width=.2) + theme_classic()+xlab("Soils")+ylab("Number of Core Taxa Observed per Sample")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))
p2

Graph relative abundance of core taxa per sample (Figure 5D)

library(Rmisc)
library(ggplot2)
Diversity_stat <- summarySE(Div, measurevar="Core_rel_abund", groupvars=c("Soils"), na.rm = TRUE)
color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
p2=ggplot() + geom_point(data=Diversity_stat, aes(x=reorder(Soils, Core_rel_abund), y=Core_rel_abund, color=Soils), size=4)+geom_point(data=Div, aes(x=Soils, y=Core_rel_abund, color=Soils), alpha=0.8) + geom_errorbar(data=Diversity_stat, aes(x=Soils, ymin=Core_rel_abund-se, ymax=Core_rel_abund+se, color=Soils), width=.2) + theme_classic()+xlab("Soils")+ylab("Relative abundance of Core Taxa ")+scale_color_manual(values=color)+ theme(axis.title = element_text(color="black", size=14, face="bold"))+ theme(axis.text = element_text(color="black", size=12, face="bold"))+ theme(legend.text = element_text(colour="black", size = 12, face = "bold"))+ theme(legend.title = element_text(colour="black", size=14, face="bold"))+ scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 0.8), breaks = c(0, 0.25, 0.5, 0.75))+ annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.4786, ymax = 0.5369, fill = "#AEBD38", alpha = .2, color = NA)+geom_hline(yintercept=0.5077, color="#AEBD38") 
p2

Figure 6

Heatmap relative abundance of core taxa (Figure 6)

library(pheatmap)
library(RColorBrewer)
library(viridis)
metacore<-read.table("Core_taxa_meta_core.txt", header=TRUE, sep = "\t")
matrix<-read.table("18S&16S_SVtable_heatmap_relabund.txt", header=TRUE, sep = "\t", row.names = 1)
meta2<-read.table("metadata_core.txt", header=TRUE)

Scale the rows

cal_z_score <- function(x){
  (x - mean(x)) / sd(x)
}
 
data_norm <- t(apply(matrix, 1, cal_z_score))

Heatmap core taxa

my_sample_col <- meta2$Soils
mat_col <- data.frame(Soils = my_sample_col)

row.names(mat_col) <- colnames(data_norm)
#faire couleur par groupe taxo (fungi...)
my_gene_col=metacore$Group2
mat_row <- data.frame(Group = my_gene_col)
row.names(mat_row) <- row.names(data_norm)
color= c("#8B5742", "#FFA07A", "#CAE2FF", "#1632AF", "#CCC0DA", "#60497A", "#C4D79B", "#6E8B3D")
ann_colors = list(
Soils = c(CAM1="#8B5742",CAM2="#FFA07A", FR1="#CAE2FF",FR2="#1632AF", IT1="#CCC0DA", IT2="#60497A", SEN1="#C4D79B", SEN2="#6E8B3D"),
Group = c( Amoebozoa="grey77", Archaea= "#1f78b4", Archaeplastida="#b2df8a", Bacteria="#1E1E1E", Hacrobia="#ff97a2",Fungi="#ff7f00", Rhizaria="#ffff00", Stramenopiles="darkmagenta", Unassigned="#cab2d6"))

pheatmap(data_norm,  annotation_col = mat_col, annotation_row = mat_row, annotation_colors = ann_colors, fontsize_row = 4,cutree_rows = 2, col=brewer.pal(9,"OrRd"), show_rownames=FALSE, show_colnames = FALSE)

Figure 7

Load packages

library(phyloseq)
library(RColorBrewer) 
library(SpiecEasi) # Network analysis for sparse compositional data  
library(network)
library(intergraph)
library(ggnet)
library(igraph)
library(microbiome)
library(ggpubr)

Import SV table and taxonomy of the core taxa (n=179 ESVs)

otu.core <- read.table(file="Core_Taxa_SVtable_for_network.txt", sep='\t', header=TRUE,check.names=FALSE,row.names=1)
taxo.core <- read.table(file="Core_Taxa_taxo_for_network.txt", sep='\t', header=TRUE,check.names=FALSE,row.names=1)
otuall3=as.matrix(t(otu.core))
otuall4=as.matrix(otu.core)
taxonomy=as.matrix(taxo.core)
TAXall=tax_table(taxonomy)
OTUall=otu_table(otuall3,taxa_are_rows=TRUE)
physeq_all = phyloseq(OTUall, TAXall)
physeq_all
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 179 taxa and 60 samples ]
## tax_table()   Taxonomy Table:    [ 179 taxa by 10 taxonomic ranks ]

Calculate network on the 179 Core taxa using the SpiecEasi package

#net.c <- spiec.easi(otuall4, method='mb', lambda.min.ratio=5e-4, nlambda=60,icov.select.params=list(rep.num=99, ncores=7))

#export
#saveRDS(net.c, "network_core16&18S_final.rds")

The output of spiec.easi is stored as network_core16&18S_final.rds

net.c <- readRDS("network_core16&18S_final.rds")
class(net.c)
## [1] "pulsar.refit"
n.c <- symBeta(getOptBeta(net.c))

Check stabilty of network (need to be close to 0.05). If not This problem we can be fixed by lowering lambda.min.ratio to explore denser networks.To get closer to the mark, we should bump up nlambda to more finely sample the lambda path, which gives a denser network.

getStability(net.c)
## [1] 0.04723226
getOptLambda(net.c)
## [1] 0.3932121
sum(getRefit(net.c))/2
## [1] 496

Prepare data for plotting using igraph package

#Add names to IDs
#We also add abundance values to vertex (=nodes=SVs).

colnames(n.c) <- rownames(n.c) <- colnames(otuall4)
# add log abundance as properties of vertex/nodes
vsize <- log2(apply(otuall4, 2, mean))

ig <- graph.adjacency(n.c, mode='undirected', add.rownames = TRUE, weighted = TRUE)
ig # we can see all the attributes and weights
## IGRAPH 7aded44 UNW- 179 496 -- 
## + attr: name (v/c), TRUE (v/c), weight (e/n)
## + edges from 7aded44 (vertex names):
## [1] 04ad657dd803395528646287365f099b--1ae16fc4b8053873d53ab4260c9311e0
## [2] 04ad657dd803395528646287365f099b--41c947a6d9a65627effb81702161b3bf
## [3] 04ad657dd803395528646287365f099b--52ef822b084fa0671f0ba6cd2b71c9ea
## [4] 04ad657dd803395528646287365f099b--6f429de2c83ae52f7700f1f348f1d88c
## [5] 04ad657dd803395528646287365f099b--ac72d8309c00939db305200bd87856c0
## [6] 04ad657dd803395528646287365f099b--ae8463313c31ed1c4fa197c805ceb8db
## [7] 04ad657dd803395528646287365f099b--bac7f39772065f7af6336d5c2beafff3
## [8] 04ad657dd803395528646287365f099b--e5ba8f3c253dc64b003e1049a2cd9de9
## + ... omitted several edges

Calculate centrality index for each node (ESV)

library(qgraph)
library(erer)
allindex=centrality_auto(ig, weighted=TRUE)

#print(allindex)

Graph relationship Node degree with Betweeness centrality or with Closeness centrality to identify hub/keystone species

index=read.table(file="Network_index_core_16S&18Sfinal.txt", header=TRUE)
taxoall3 <- read.table(file="Core_Taxa_taxo_for_network.txt", sep='\t', header=TRUE,check.names=FALSE)

index_tax=merge(index, taxoall3, by="OTUID")

Graph Node degree vs Betweeness centrality (Figure 7C)

library(ggplot2)
mycolors <- scale_color_manual(values = c("#a6cee3",  "#b2df8a","#1f78b4","#1E1E1E","#000075","#33a02c","#fdbf6f","#ff7f00","#cab2d6","#ffff99","#b15928",'#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#000000',"#800000","#8E388E","#7171C6","#7D9EC0","#388E8E","#71C671","#8E8E38","#C5C1AA", "#C67171","#555555", "orange"))
myshape <- scale_shape_manual(values=c(17, 16))
p1=ggplot(data=index_tax, aes(x=Degree, y=Betweenness,color=Class, shape=Hub.x))+geom_point(size=3)+theme_classic(base_size = 8)+xlab("Node Degree (nb Correlations)")+ylab("Betweenness Centrality ")+ theme(axis.title = element_text(color="black", size=12, face="bold"))+ theme(axis.text = element_text(color="black", size=10, face="bold"))+ theme(legend.text = element_text(colour="black", size = 6, face = "bold"))+ theme(legend.title = element_text(colour="black", size=8, face="bold"))+mycolors+myshape
p1

Graph Node degree vs Closeness centrality (Figure 7D)

library(ggplot2)
mycolors <- scale_color_manual(values = c("#a6cee3",  "#b2df8a","#1f78b4","#1E1E1E","#000075","#33a02c","#fdbf6f","#ff7f00","#cab2d6","#ffff99","#b15928",'#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#000000',"#800000","#8E388E","#7171C6","#7D9EC0","#388E8E","#71C671","#8E8E38","#C5C1AA", "#C67171","#555555", "orange"))

p1=ggplot(data=index_tax, aes(x=Degree, y=Closeness,color=Class, shape=Hub.x))+geom_point(size=3)+theme_classic(base_size = 8)+xlab("Node Degree (nb Correlations)")+ylab("Closeness Centrality ")+ theme(axis.title = element_text(color="black", size=12, face="bold"))+ theme(axis.text = element_text(color="black", size=10, face="bold"))+ theme(legend.text = element_text(colour="black", size = 6, face = "bold"))+ theme(legend.title = element_text(colour="black", size=8, face="bold"))+mycolors+myshape
p1

Improved network visualization with ggnet package, add colors to edges

net <- asNetwork(ig)
network::set.edge.attribute(net, "color", ifelse(net %e% "weight" > 0, "lightblue", "red"))

add information on taxonomy, hub, nodesize for each core taxon in the network

class <- map_levels(colnames(otuall4), from = "OTUID2", to = "Supergroup", tax_table(physeq_all))
net %v% "Supergroup" <- class
cluster <- map_levels(colnames(otuall4), from = "OTUID2", to = "Cluster", tax_table(physeq_all))
net %v% "Cluster" <- cluster
net %v% "nodesize" <- vsize

hub2 <- index_tax$Hub.x
hub2 <-  as.character(hub2)


net %v% "Hub" <- hub2

Network graph - nodes colored by Taxonomic Group (Figure 7A)

mycolors <- scale_color_manual(values = c("grey77",  "#1f78b4","#b2df8a","#1E1E1E","#C67171","#ff7f00","pink","goldenrod1","darkmagenta","#cab2d6","#ffff99","#b15928",'#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#000000',"#800000","#8E388E","#7171C6","#7D9EC0","#388E8E","#71C671","#8E8E38","#C5C1AA", "#C67171","#555555", "orange"))
myshape <- scale_shape_manual(values=c(17, 16))

p <- ggnet2(net, mode = "fruchtermanreingold", layout.par = list(cell.jitter = 0.75), node.color = "Supergroup", label = FALSE, node.size = "nodesize", edge.color = "color", node.shape = "Hub", max_size = 4, legend.size = 12) + guides(color=guide_legend(title="Group"), size = FALSE, mode = c("x", "y")) + mycolors +myshape
p 

Network graph - nodes colored by core taxa cluster (Figure 7B)

mycolors <- scale_color_manual(values = c("darkblue",  "darkred","gray"))
myshape <- scale_shape_manual(values=c(17, 16))

p <- ggnet2(net, mode = "fruchtermanreingold", layout.par = list(cell.jitter = 0.75), node.color = "Cluster", label = FALSE, node.size = "nodesize", edge.color = "color", node.shape = "Hub", max_size = 4, legend.size = 12) + guides(color=guide_legend(title="Cluster"), size = FALSE, mode = c("x", "y")) + mycolors+myshape
p